//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************

#include "PrimaryGeneratorAction.hh"

#include <G4Types.hh>

#include "DetectorConstruction.hh"
#include "G4ParticleDefinition.hh"
#include "G4ParticleGun.hh"
#include "G4SystemOfUnits.hh"
#include "G4Threading.hh"
#include "MyMacros.hh"
#include "PrimaryGeneratorMessenger.hh"
#include "Randomize.hh"

// #include "G4ParticleTable.hh"
// #include "G4PhysicalConstants.hh"

constexpr G4double _beamEnergyDft{480.0};

PrimaryGeneratorAction::PrimaryGeneratorAction(DetectorConstruction* dtctr) : _fDetector(dtctr)
{
#if 0
    MYG4TID(tid);
    DINFOLN("```````````````````````````````````````````[tid {} PrimaryGeneratorAction/PrimaryGeneratorAction]", tid);
#endif
    InitializeMe();
}

void PrimaryGeneratorAction::InitializeMe()
{
    _fVerbose = _fDetector->GetVerbose();

    _fMessenger = new PrimaryGeneratorMessenger(this);
    _fParticleGun = new G4ParticleGun();

    _fPosition = G4ThreeVector(_fX0, _fY0, _fZ0);
    _fDirection = G4ThreeVector(0.0, 0.0, 1.0);

    SetBeamEnergy(_beamEnergyDft * MeV);
}

void PrimaryGeneratorAction::ReInitialize()
{
    _fPosition = G4ThreeVector(_fX0, _fY0, _fZ0);
    _fDirection = G4ThreeVector(0.0, 0.0, 1.0);
    SetBeamEnergy(_beamEnergyDft * MeV);
    _fCounter = 0;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

PrimaryGeneratorAction::~PrimaryGeneratorAction()
{
    delete _fParticleGun;
    delete _fMessenger;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
#if 0
    MYG4TID(tid);
    DINFOLN("[tid {} PrimaryGeneratorAction] GeneratePrimaries()", tid);
#endif
    _fCounter++;

    // Simulation of beam position
    G4double x = _fX0;
    G4double y = _fY0;
    G4double z = _fZ0;  //2022.10.18
    _fPosition = G4ThreeVector(x, y, z);
    _fParticleGun->SetParticlePosition(_fPosition);

    // for QA beam model
    G4double spoty = _fSpotDeltaY;
    G4ThreeVector posVSCenter = G4ThreeVector(0.0, spoty, 0.0);  // position of keypoint (between source and iso-center, in the spot-plane)

    _fRMax2 = std::pow(5.0 * std::max(_fSigmaX, _fSigmaY), 2);  //2022.10.14
    G4double dx = 0, dy = 0;
    if (_fDisTypeXY == "flat")
    {
        dx = _fSigmaX * (1 - 2 * G4UniformRand());
        dy = _fSigmaY * (1 - 2 * G4UniformRand());
    }
    else if (_fDisTypeXY == "gauss")
    {
        do
        {
            if (0.0 < _fSigmaX)
            {
                dx = G4RandGauss::shoot(0.0, _fSigmaX);
            }
            if (0.0 < _fSigmaY)
            {
                dy = G4RandGauss::shoot(0.0, _fSigmaY);
            }
        } while (dx * dx + dy * dy > _fRMax2);
    }

    G4ThreeVector locPos = G4ThreeVector(dx, dy, 0.);
    G4ThreeVector dirX = G4ThreeVector(1.0, 0.0, 0.0);
    G4ThreeVector dirZ = (posVSCenter - _fPosition).unit();
    G4ThreeVector dirY = dirZ.cross(dirX);
    G4ThreeVector posVS = posVSCenter + RotationMatrix(locPos, dirX, dirY, dirZ);

    _fDirection = (posVS - _fPosition).unit();

#if 0
  // Simulation of beam direction
  G4double ux = fDirection.x();
  G4double uy = fDirection.y();
  G4double uz = fDirection.z();

  // Beam particles are uniformly distributed over phi, cosTheta
  if(fDisTypeTheta == "flat") {
      G4double minCosTheta = std::cos(fSigmaTheta);
      uz = minCosTheta + (1.0 - minCosTheta)*G4UniformRand() ;
      ux = std::sqrt((1.0 - uz)*(1.0 + uz)) ;
  } else if(fDisTypeTheta == "gauss") {
      if(0.0 < fSigmaTheta){
          uz = std::cos(G4RandGauss::shoot(0.0,fSigmaTheta));
          ux = std::sqrt((1.0 - uz)*(1.0 + uz));
      }
  }

  G4double phi = twopi*G4UniformRand() ;
  uy = ux;
  ux *= std::cos(phi) ;
  uy *= std::sin(phi) ;
  fDirection = G4ThreeVector(ux,uy,uz) ;
#endif

    _fParticleGun->SetParticleMomentumDirection(_fDirection);

    ///======================= Simulation of beam kinetic energy
    G4double kinEnergy = _fEnergy;
    if (_fDisTypeE == "flat")
    {
        kinEnergy = _fEnergy + _fSigmaE * _fEnergy * (1 - 2 * G4UniformRand());
    }
    else if (_fDisTypeE == "gauss")
    {
        if (0.0 < _fSigmaE) kinEnergy = _fEnergy + G4RandGauss::shoot(0.0, _fSigmaE * _fEnergy);
    }

    _fParticleGun->SetParticleEnergy(kinEnergy);

    if (G4Threading::G4GetThreadId() == 0 && _fCounter == 1)
    {
        PrintSourceParameters();
    }

    if (_fVerbose > 0)
    {
        G4ParticleDefinition* particle = _fParticleGun->GetParticleDefinition();
        G4String particleName = particle->GetParticleName();
#if 0
        G4cout << "Event#  " << _fCounter << "  Beam particle is generated by PrimaryGeneratorAction " << G4endl;
        G4cout << "ParticleName= " << particleName << "  PDGcode= " << particle->GetPDGEncoding() << std::setprecision(5)
               << "   KinEnergy(GeV)= " << _fEnergy / MeV << "   x(cm)= " << x / cm << " y(cm)= " << y / cm << " z(cm)= " << z / cm << G4endl;

#endif 0
        LINFO(
            "```````````````````````````"
            "[PrimaryGeneratorAction/GeneratePrimaries] "
            "Event {}, Particle {}, PDGcode {}, KinEngy(MeV) {:>11.6F}, "
            "x(cm) {}, y(cm) {}, z(cm) {}",
            _fCounter, particleName, _fEnergy / MeV, x / cm, y / cm, z / cm);
    }

    _fParticleGun->GeneratePrimaryVertex(anEvent);
}

void PrimaryGeneratorAction::SetBeamEnergy(G4double val) { _fEnergy = val; }

void PrimaryGeneratorAction::PrintSourceParameters()
{
    MYG4TID(tid);
    G4ParticleDefinition* particle = _fParticleGun->GetParticleDefinition();
    G4String particleName = particle->GetParticleName();

    // G4cout << "===================================================" << G4endl;
    // G4cout << "#              Geant4 test Source                 #" << G4endl;
    // G4cout << "===================================================" << G4endl;
    LINFO("``````````````````````````````````````````````````[tid {} PrimaryGeneratorAction/PrintSourceParameters]", tid);
    G4cout << "  Event#  " << _fCounter << G4endl;
    G4cout << "  particle name= " << particleName << G4endl;
    G4cout << "  PDG code= " << particle->GetPDGEncoding() << G4endl;
    G4cout << "  central energy= " << _fEnergy / MeV << " MeV " << G4endl;
    G4cout << "  energy dis. type= " << _fDisTypeE << ", sigma(rel.)= " << _fSigmaE << G4endl;
    G4cout << "  central position= (" << _fX0 / mm << ", " << _fY0 / mm << ", " << _fZ0 / mm << ") mm " << G4endl;
    G4cout << "  pos. dis. type= " << _fDisTypeXY << ", sigma= (" << _fSigmaX / mm << ", " << _fSigmaY / mm << ", " << _fSigmaZ / mm
           << ") mm, RMax= " << std::sqrt(_fRMax2) / mm << " mm " << G4endl;
    G4cout << "  initial direction= (" << _fDirection.x() << ", " << _fDirection.y() << ", " << _fDirection.z() << ") " << G4endl;
    G4cout << "  central theta=" << 0. << G4endl;
    G4cout << "  theta dis. type= " << _fDisTypeTheta << ", sigma= " << _fSigmaTheta / rad << " rad " << G4endl;
    G4cout << "===================================================" << G4endl;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

inline G4ThreeVector PrimaryGeneratorAction::RotationMatrix(G4ThreeVector& locPos,
                                                            G4ThreeVector& Rotx,
                                                            G4ThreeVector& Roty,
                                                            G4ThreeVector& Rotz) const
{  //2023.6.20
    G4double x = locPos.x() * Rotx.x() + locPos.y() * Roty.x() + locPos.z() * Rotz.x();
    G4double y = locPos.x() * Rotx.y() + locPos.y() * Roty.y() + locPos.z() * Rotz.y();
    G4double z = locPos.x() * Rotx.z() + locPos.y() * Roty.z() + locPos.z() * Rotz.z();
    return G4ThreeVector(x, y, z);
}
